#!/usr/bin/env python
"""
Reads an HMM file and converts each natural log probability into a position specific
scoring matrix.

"""
import math
import re
import os
import sys
import multiprocessing
from Bio import SeqIO, AlignIO
import time
import subprocess
#############    < OPTIONS >	################

options = { 	 '-hmm' : [] ,
		 '-ncpus' : multiprocessing.cpu_count()-1,
		 '-out' : sys.stdout,
		 '-in'  : None,
		 '-format':'fasta',
		 '-desc': [],
		 '-accs': [],
		 '-hmmalign' : '/biol/bin/hmmalign',
		 '--consensus':False
	}
##		HELP			#####
helpTxt = {	'-hmm' : 'HMM input file.',
		'-ncpus':'Number of processes to intitiate which will do the math calculations',
		'-out' : 'Output file. Default: standard out. N.B. Appends to the file.',
		'-in'  : 'Input sequence file. Required.',
		'-desc': 'Filter input sequences by description(s), space delimited (Enclose with \' or " ).',
		'-accs': 'Filter input sequences by accession(s), space delimited.',
		'--consensus': "Also compare each HMM consensus sequence to the model whence it came. [False]",
		'-hmmalign':'hmmalign program. Case sensitivity matters with the alignments!'
	  }


def printHelp():
	global helpTxt
	print 'Valid options:-'
	for key,val in helpTxt.items():
		print key,' : ',val

def parseArgs(args):
	global options
	global helpTxt
	if len(args) == 0:
		printHelp()
		exit()
	prevargDash = False
	for i in xrange(len(args)):
		if args[i].startswith('--'):
			if args[i] in options.keys():
				options[ args[i] ] = True
			else:
				print 'Invalid option: {0}'.format( args[i] )
				printHelp()
				exit()
		elif args[i].startswith('-'):
			if args[i] in options.keys():
				prevargDash = True
				command = args[i]
			else:
				print 'Invalid option: {0}'.format( args[i] )
				printHelp()
				exit()
		elif prevargDash:
			if command in ['-in','-out']:
				options.update( { command : args[i] } )
				prevargDash = False
			elif command == '-accs':
				options[command].append(args[i])
			elif command == '-desc':
				options[command].append( '('+args[i] + ')')
			elif command == '-format':
				options.update( {command : args[i] } )
				prevargDash = False
			elif command == '-hmm':
				options[command].append( args[i] )
		else:
			print "Don't know what to do with {0}".format( args[i] )
			exit()
	if options['-in'] == None:
		sys.stderr.write( '"-in" is a required argument!\n')
		exit()
	return options

#############    </HELP & OPTIONS>	################

class ColError( Exception ):
	def __init__(self):
		Exception.__init__(self,"\nAlignment column order wrong in HMM!!")


############	MEAT & VEG	##############
def parseHMM( options, prob_calc_procs ):
	if options['-hmm'] == None:
		print 'No HMM specified'
		for proc in prob_calc_procs:
			proc.inQ.put( 'endSEQ' )
			proc.inQ.put( 'STOP' )
			proc.join()
		exit()
	else:
		inHandle = file(options['-hmm'],'r')
	split_cols = re.compile( r'\s+' )
	line = inHandle.readline()
	HMMLineReg = re.compile(r'HMM\s+') ## Finds line that starts 'HMM' followed by white space.
	while not HMMLineReg.search(line): ## Go to that line. (Contains column headers)
		line = inHandle.readline()
	i = 0
	matchCols = {}
	for residue in split_cols.split( line.rstrip() )[1:]:  ## Column headers. For residue in header line:
		matchCols.update( { residue : i } )	       #### Add { residue : column# } to matchCols.
		i += 1
	inHandle.readline()	    ## Line underneath HMM line. Ignore.
	globalProbs = split_cols.split( inHandle.readline().strip() )[1:]   ## The COMPO line = Overall probabilities per residue.
	insertProbs = [split_cols.split( inHandle.readline().strip() )]	    ## Insert probabilities at position 0 (before HMM).
				    #7 entries. Transition from begin / insert state 0
	begins = split_cols.split( inHandle.readline().strip() )
				    # B-->M1 ; B-->I0 ; B-->D1 ; I0-->M1 ; I0-->I0; P=1 ; P=0.
				    # ..Must use one.
	line = inHandle.readline()	# First proper line.
	con_col = 1  # Counter for consensus columns. HMMer asserts column order with first entry in each group of 3 lines.
	matchProbs = [ ]
	transitionProbs = [ ]
	while line != '//':
		entries = split_cols.split( line.strip() )
		try:
			assert int(entries[0]) == con_col
		except AssertionError:
			raise ColError(entries[0])
		matchProbs.append( entries[1:] ) ## Making the array. Hiccup: indexes are 1 No. below column No. 1 per residue & Col# ; RF Annotation ; Consensus structure.
		insertProbs.append( split_cols.split( inHandle.readline().strip() ) ) # 1 per residue.
		transitionProbs.append( split_cols.split( inHandle.readline().strip() ) ) # Mn-->Mn+1 ; In ; Dn+1 ; In-->Mn+1 ; In ; Dn-->Mn+1 ; Dn+1
		try:
			assert len( transitionProbs[-1] ) == 7
		except Exception:
			sys.stderr.write( 'Transition line (below) does not have 7 states...' )
			sys.stderr.write( transitionProbs[-1] )
			raise
		line = inHandle.readline().strip()
		con_col += 1
	matrices = {	'col_index' : matchCols,
			'match' : matchProbs,
			'insert': insertProbs,
			'globals': globalProbs,
			'begin'	: begins,
			'transitions' : transitionProbs
			}
	inHandle.close()
	return matrices

class calc_prob( multiprocessing.Process ):
	def __init__(self,inQ,sumQ,lock):
		multiprocessing.Process.__init__(self)
		self.inQ = inQ
		self.sumQ = sumQ
		self.lock = lock
	def run( self ):
		in_LO = self.inQ.get() ## Log Odds ratio parsed from HMM.
		while in_LO != 'STOP':
			while in_LO != 'endSEQ' and in_LO != 'STOP':
				try:
					res, ResScores = in_LO.items()[0]
					NSeqs = len( ResScores)
					if NSeqs == 1:
						Prob = math.exp( - float(ResScores[0]) ) 
					elif NSeqs > 1:
						Prob = math.fsum([math.exp( - float(p) ) for p in ResScores ] ) / float( NSeqs ) ## Average probabilities
					LogOdds = - math.log( 1 - Prob ) ## Log Odds probability of residue not being observed.
					self.sumQ.put( {res : LogOdds } )
					in_LO = self.inQ.get()
				except Exception:
					sys.stderr.write( str(in_LO ))
					sys.stderr.flush()
					self.sumQ.put( 'endSEQ' )
					self.sumQ.put('STOP' )
					raise
			if in_LO == 'endSEQ':
				time.sleep(0.01)
				self.sumQ.put( in_LO ) ## Send 'endSEQ'
				in_LO = self.inQ.get()

class sum_probs( multiprocessing.Process ):
	def __init__(self,inQ,outQ ):
		multiprocessing.Process.__init__(self)
		self.inQ = inQ
		self.outQ = outQ
	def run(self):
		in_value = self.inQ.get()
		while in_value != 'STOP':
			tally = float(0)
			probs = {}
			while in_value != 'endSEQ' and in_value != 'STOP': 
				res, LogOdds = in_value.items()[0]
				tally += LogOdds   ## Sum the Log Odds probabilities.
				if res in probs.keys():
					probs[res].append(LogOdds)
				else:
					probs.update( {res : [LogOdds]} )
				in_value = self.inQ.get()
			if in_value == 'STOP':
				break
			else:
				self.outQ.put( [tally,probs] )
				in_value = self.inQ.get()
		self.inQ.close()
		self.outQ.close()

def score_alignment( alignment_residues, matrices, processes ):
	IUPACcodes = { 	##'T' : ['U'], ## Will consistently only turn U --> T.
			'T' : ['T'],
			'A' : ['A'],
			'G' : ['G'],
			'C' : ['C'],
			'U' : ['T'],
			'M' : ['A','C'],
			'R' : ['A','G'],
			'W' : ['A','T'],
			'S' : ['C','G'],
			'Y' : ['C','T'],
			'K' : ['G','T'],
			'V' : ['A','C','G'],
			'H' : ['A','C','T'],
			'D' : ['A','G','T'],
			'B' : ['C','G','T'],
			'X' : ['G','A','T','C'],
			'N' : ['G','A','T','C']
			}
	col_ind = 0
	proc_end = len(processes)
	proc_ind = 0
	residues = matrices['col_index']
	if 'U' in residues:
		residues.update( { 'T' : residues['U'] } )
	state = 'begin'
	col_ind = 0
	res_ind = 0
	seq = alignment_residues
	## B --> M1, B --> I0, B --> D1; I0 --> M1, I0 --> I0; then a 0.0 and a .*.
	while state != 'match' and state != 'delete':  ## Look through until leave begin state.
		probs = []
		res = seq[res_ind]
		if res in ['-','~']:		## Delete state.
			processes[proc_ind].inQ.put( { col_ind :[matrices['begin'][2]] } )
			proc_ind += 1
			state = 'delete'
			col_ind += 1
		elif res in ['.']:	## Insert state
			if state == 'insert':
				processes[proc_ind].inQ.put( { col_ind : matrices['begin'][4] } )
			elif state == 'begin':
				processes[proc_ind].inQ.put( { col_ind : matrices['begin'][1] } )
			state = 'insert'
			proc_ind += 1
		elif res.isupper():	## Match state.
			assert col_ind == 0
			if state == 'begin':
				processes[proc_ind].inQ.put( {col_ind :[matrices['begin'][0]] } )  ## B --> M1
			elif state == 'insert':
				processes[proc_ind].inQ.put( {col_ind :[matrices['begin'][3]]} )  ## I0 --> M1
			else:
				raise ValueError("The state {0} is not valid for before the first match state".format(state))
			for trans in IUPACcodes[res.upper() ]:
				probs.append( matrices['match'][col_ind][ residues[ trans] ] )
			processes[proc_ind].inQ.put( { col_ind : probs } )
			state = 'match'
			col_ind += 1
			proc_ind += 1
		else:			## Lower case == insert state.
			if state == 'insert':
				processes[proc_ind].inQ.put( { col_ind : [ matrices['begin'][4] ] }) ## I0 --> I0
			elif state == 'begin':
				processes[proc_ind].inQ.put( { col_ind : [ matrices['begin'][1] ] }) 	## B-->I0
			else:
				sys.stderr.write( "Can't go from state {0} to insert".format(state) )
			state = 'insert'
			for trans in IUPACcodes[res.upper() ]:
				probs.append( matrices['match'][col_ind][ residues[trans] ] )
			processes[proc_ind].inQ.put( {col_ind :probs} )
		if proc_ind >= proc_end :
			proc_ind = 0
		res_ind += 1

	for res in seq[col_ind:]:
		try:
			probs = []
			if res in [ '-', '~'  ]: # Allowed gap residues (and probably more...).
				if state == 'match':
					processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][2]] } ) ## Mk --> Dk+1
				elif state == 'delete':
					delVal = matrices['transitions'][col_ind][6]
					if delVal != '*':
						processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][6]] } ) ## Dk --> Dk+1
				else:
					sys.stderr.write( "Can't go from insert state to delete state!" )
					sys.stderr.flush()
				state = 'delete'
				col_ind += 1
			elif res == '.':  	# Insert state
				if state == 'insert':
					in_prob = matrices['transitions'][col_ind][4]
				elif state == 'match':
					in_prob = matrices['transitions'][col_ind][1]
				else:
					sys.stderr.write( "Can't go from state {0} to 'insert'".format(state) )
				processes[proc_ind].inQ.put( {col_ind : [in_prob] } )
				for trans in IUPACcodes[res.upper() ]:
					probs.append( matrices['insert'][col_ind][ residues[trans] ] )
				state = 'insert'
			elif res.isupper():
				if state == 'match':
					processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][0]]} ) ## Mk --> Mk+1
				elif state == 'delete':
					processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][5]]} ) ## Dk --> Mk+1
				elif state == 'insert':
					processes[proc_ind].inQ.put( {col_ind : [matrices['transitions'][col_ind][3]]} ) ## Ik --> Mk+1
				else:
					sys.stderr.write( "Unknown state: {0}.\n".format(state) )
				for trans in IUPACcodes[ res ]:
					probs.append( matrices['match'][col_ind][ residues[trans ] ])
				state = 'match'
				col_ind += 1
			else:
				if state == 'insert':
					in_prob = matrices['transitions'][col_ind][4]  ## Ik --> Ik
				elif state == 'match':
					in_prob =  matrices['transitions'][col_ind][1]  ## Mk --> Ik
				else:
					sys.stderr.write( "Can't go from state {0} to 'insert'".format(state) )
				processes[proc_ind].inQ.put( {col_ind : [in_prob] } )
				try:
					for trans in IUPACcodes[res.upper() ]:
						probs.append( matrices['insert'][col_ind][ residues[trans] ] )
				except Exception:
					print col_ind
					print len( matrices['insert'] )
					print residues
					print matrices['insert'][col_ind]
					raise
				state = 'insert'
			proc_ind += 1
			if proc_ind >= proc_end:
				proc_ind = 0
			processes[proc_ind].inQ.put( {col_ind :probs } ) ## -ve log odds probability.
			res_ind += 1
		except Exception:
			print '\ncol ind',col_ind
			print 'len insert: ',len( matrices['insert'] )
			print 'len matches: ',len( matrices['match'] )
			print 'len transi: ',len(matrices['transitions'])
			print 'seqlen: ',len(seq)
			print seq
			print 'residues',residues
			print 'matrices[insert][col_ind]:,', matrices['insert'][col_ind]
			#print 'matrices[match][col_ind]:,', matrices['match'][col_ind]
			print 'proc_ind:',proc_ind
			print 'proc_end:',proc_end
			raise
	processes[proc_ind].inQ.put( 'endSEQ' )
	return True

#############	</MEAT & VEG>	############

def get_consensus( HMM ):
	emit_process = subprocess.Popen( [ 'hmmemit','-c', HMM ], shell = False, stdout = subprocess.PIPE )
	consensus_seq = SeqIO.read( emit_process.stdout , 'fasta' )
	emit_process.stdout.close()
	return consensus_seq


def parse_seqs(options):
	if type(options['-in']) != file:
		options.update( { '-in' : file(options['-in'],'r') } )
	if len(options['-desc']) > 0:
		desc_reg = re.compile( '|'.join([arg for arg in options['-desc']] ), re.I )
	noFilter = False
	if len(options['-accs']) > 0:
		acc_reg = re.compile( '|'.join([re.escape(arg) for arg in  options['-accs']] ), re.I )
	elif len(options['-accs']) == 0 and len(options['-desc']) == 0:
		noFilter = True
	Ndesc = len(options['-desc'])
	Naccs = len(options['-accs'])
	count = 0
	for seq in SeqIO.parse( options['-in'],options['-format']  ):
		count += 1
		if Naccs > 0:
			accFind = acc_reg.search( seq.id )
			if accFind:
				sys.stdout.flush()
				yield seq
				continue
		if Ndesc > 0:
			if desc_reg.search( seq.description ):
				yield seq
		elif noFilter:
			yield seq

class align_seq( multiprocessing.Process ):
	def __init__(self,inPipe,outPipe,options):
		multiprocessing.Process.__init__(self)
		self.in_pipe = inPipe
		self.out_pipe = outPipe
		self.options = options.copy()
		self.options.update( { '-acc':[] , '-desc':[], '-informat':options['-format']  } )
		self.options.update( { '-format' : 'stockholm' } )
		#self.sem = semaphore
		self.hmmalign = options['-hmmalign']
		self.hmmsearch = 'hmmsearch'
	def run(self):
		inp = self.in_pipe.recv()
		while inp != 'STOP':
			seq_record , HMM = inp
#			align_process = subprocess.Popen( [ self.hmmalign,'--informat',self.options['-informat'],'--outformat','A2M',HMM,'/dev/stdin'],shell=False,stdin=subprocess.PIPE,stdout=subprocess.PIPE )
			align_process = subprocess.Popen( [ self.hmmsearch,'-A','/dev/stdout','-o',os.devnull,HMM,'/dev/stdin'],shell=False,stdin=subprocess.PIPE,stdout=subprocess.PIPE )
			align_process.stdin.write( '>{0}\n'.format( seq_record.id ))
			align_process.stdin.write( seq_record.seq.tostring().replace( '-','' ) )
			options.update( { '-in' : align_process.stdout } )
			align_process.stdin.close()
			for seq in AlignIO.read( align_process.stdout , 'stockholm' ):
				self.out_pipe.send_bytes( seq.seq.tostring() )
			align_process.wait()
			align_process.stdout.close()
			inp = self.in_pipe.recv()
		self.in_pipe.close()
		self.out_pipe.close()
		
def print_seqs_vs_HMM( accession_names , matrices , out = sys.stdout ):
	"""This prints multiple seqeunces against each sequential HMM. So the output
	will contain columns containing the cumulative sequence scores of each sequence
	against just one model. If there's multiple HMMs, then the next HMM will be
	printed underneath.
	"""
	HMMs = sorted( matrices.keys() )
	accessions = sorted(accession_names.keys())
	out.write('\n')
	out.write('# Each section is split by headers where each line represents a different descriptive header.\n')
	out.write('# The first line gives the name of the HMM.\n')
	out.write('# The second line always the accession numbers for the target sequences.\n')
	out.write('# The next line will give the first 30 characters in the descriptive sequence header.\n')
	out.write('# After that comes all the marvelous data....\n')
	out.write('\n# These HMMs will be searched, in this order: ')
	out.write( str(', '* len(accessions)).join(  HMMs ) +'\n\n')
	if 'cons' in accessions:
		del(accessions[ accessions.index('cons') ] )
		accessions.append('cons')
	for n in xrange( len(HMMs) ):
		## Print column headers
		out.write( '# ***|' + str('\t'*len(accessions)) + HMMs[n] + str('\t'*len(accessions))+'|***\n')
		out.write( '#\t'+ '\t'.join( [accession_names[acc][:30] for acc in accessions ] ) +'\n')  # Sequence description line
		out.write('#Col No.\t'+ ( '\t'.join(accessions))+'\n')					## Sequence accession
		tallies = [ float(0) for i in xrange( len(accessions ) )]
		matrix = matrices[HMMs[n]]
		print [len(matrix[acc][1].keys()) for acc in accessions]
		lengths = [ len(matrix[acc][1].keys()) for acc in accessions ]
		nResidues = max( lengths )
		for i in xrange( nResidues):
			out.write('\n{0}\t'.format(i))
			colInds = {}
			inserts = []
			for seqInd in xrange(len(accessions)):
				if lengths[seqInd] <= i:
					matrix[ accessions[seqInd] ][1].update( { i: [] } )
				for score in matrix[accessions[seqInd]][1][i]:
					tallies[seqInd]  += score
			out.write( '\t'.join([str(t) for t in tallies] ) )
		out.write('\n\n')

def print_HMMs_vs_seq( accession_names , matrices , resultseq, out=sys.stdout ):
	"""This will write (to out) the cumulative sequence scores a single sequence against
	multiple HMMs. The columns are subdivided into groups, each group corresponding to an
	HMM. Within each group of columns, every sequence gets a column of it's own.
	Note that consensus sequence scores show the cumulative scores for consensus sequences
	that are unique for each HMM."""
	HMMs = sorted( matrices.keys())
	accessions = sorted(accession_names.keys())
	if 'cons' in accessions:
		del(accessions[ accessions.index('cons') ] )
		accessions.append('cons')
	out.write('# Each section is split by headers where each line represents a different descriptive header.\n')
	out.write('# The first line gives the name of the HMM.\n')
	out.write('# The second line always the accession numbers for the target sequences.\n')
	out.write('# The next line will give the first 30 characters in the descriptive sequence header.\n')
	out.write('# After that comes all the marvelous data....\n')
	out.write('\n# These HMMs will be searched, in this order: ')
	out.write( ', '.join( HMMs ) +'\n')  # Sequence description line
	out.write( '\t'.join(accessions) +'\n\n')
	if 'cons' in accessions:
		del(accessions[ accessions.index('cons') ] )
		accessions.append('cons')
	nSeqs = len(accessions)
	HMMLens = []
	maxHMMLen = 0
	for HMM in HMMs:
		HMMLens.append(  [len(matrices[HMM][accession][1].keys()) for accession in accessions]  )
		maxHMMLen = max( [maxHMMLen] + HMMLens[-1] )
	#HMMLens = [ len(matrices[HMM][accessions[0]][1].keys()) for HMM in HMMs]
	nHMMs = len(HMMs)
	ResInd = 0
	tallies = [ [float(0) for i in xrange(nSeqs)] for HMM in xrange(len(HMMs))]
	tallyInd = 0
	insertReg = re.compile( r'[a-z\.]{1}' )
	try:
		while ResInd < maxHMMLen:
			line = []
			HMMInd = 0
			while HMMInd < nHMMs:
				SeqInd = 0
				while SeqInd < nSeqs: # Each sequence gets a column in each row.
					if ResInd < HMMLens[HMMInd][SeqInd]:
						scores = matrices[ HMMs[HMMInd] ][ accessions[SeqInd] ][1][ ResInd ]
						tallies[HMMInd][SeqInd] += sum( scores )
						line.append( tallies[HMMInd][SeqInd] )
					else:
						line.append( tallies[HMMInd][SeqInd] )
					SeqInd += 1
					tallyInd = SeqInd
				HMMInd += 1
			out.write('{0}\t'.format(ResInd)+ '\t'.join( str(score) for score in line ) +'\n')
#			if insertReg.search( resultseq[ResInd] ):
#				continue
#			else:
			ResInd += 1
	except Exception:
		print HMMLens, HMMs,HMMInd, SeqInd,nSeqs,accession_names.keys(),matrices.keys()
		print accessions[SeqInd]
		print "ResInd: ",ResInd
		print "MaxHMM Len: ",HMMLens[HMMInd]
		raise
		count = 0
		for i in sorted(matrices[HMMs[HMMInd]][accessions[SeqInd] ][1]):
			try:
				assert i == count
			except AssertionError:
				sys.stderr.write( 'i: {0}, count: {1}, len: {2}\n'.format( i, count, len(matrices[HMMs[HMMInd]][accessions[SeqInd]][1])))
				raise
			count += 1
		print "N Consensus Cols: ",count
		raise



if __name__ == '__main__':
	args = sys.argv[1:]
	options = parseArgs( args )
	#### PROCESS INITIATION ####
	## Start Qs & Pipes between sum_probs & calc_prob ####
	calc_to_sumQ = multiprocessing.Queue()
#	sum_end,calc_end = multiprocessing.Pipe()
	sum_seq_queue = multiprocessing.Queue()
	## Initiate the sum_probs process ##
	sum_thread = sum_probs( calc_to_sumQ, sum_seq_queue )
	sum_thread.start()
	lock = multiprocessing.RLock()
	## Initiate the calc_prob processes ##
	math_processes = [ calc_prob(multiprocessing.Queue(),calc_to_sumQ,lock) for i in xrange(options['-ncpus']) ]
	## hmmalign controller processes & semantics ##
	seq_send_pipe, seq_recv_pipe = multiprocessing.Pipe()
	ali_recv_pipe, ali_send_pipe = multiprocessing.Pipe()
#	semaphore = multiprocessing.Semaphore(options['-ncpus'])
#	align_processes = [ align_seq( seq_recv_pipe,ali_send_pipe,semaphore,options) for i in xrange(options['-ncpus'] ) ]
	align_process =  align_seq( seq_recv_pipe,ali_send_pipe,options) 
	for i in xrange(options['-ncpus']):  ## Start the bad boiiis
		math_processes[i].start()
	align_process.start()
	####  END PROCESS INITIATION #####
	sequence_scores  = {}
	acc_names = {}
	sys.stdout.write( '# '+ '\t'.join( ['Accession', 'HMM', 'Score'] ) )
	if options['--consensus']:
		acc_names.update( { 'cons' : 'HMM Consensus sequence' } )
		sys.stdout.write( '\t{0}'.format( 'Cons. Score' ) )
	env = options.copy()
	for seq in parse_seqs( options ):  ## Iterate through sequences
		acc_names.update( { seq.id : seq.description[len(seq.id):].strip() } )
		for HMM in options['-hmm']:	## Iterate through HMMs.
			env.update( { '-hmm' : HMM } )
			matrices = parseHMM( env,math_processes )  ## Convert the HMM to log odds probability matrices.
			seq_send_pipe.send([seq,HMM])		   ## Send sequence to be aligned.
			ali = ali_recv_pipe.recv_bytes()	   ## Receive aligned seqeunce.
			submitted = score_alignment( ali, matrices, math_processes )  ## Score sequence against HMM.
			tally, probs = sum_thread.outQ.get()	   ## Get the results.
			HMMFname = HMM[HMM.rfind(os.path.sep)+1:]  ## Don't need the full path name; just the file name when writing out.
			if HMMFname in sequence_scores.keys():
				sequence_scores[HMMFname].update( { seq.id : [tally, probs] } )
			else:
				sequence_scores.update( { HMMFname : { seq.id : [tally,probs] } } )  ## Update dictionary of results.
			sys.stdout.write('\n# {0}\t{1}\t{2}'.format( seq.id,HMMFname,tally ) )
			if options['--consensus']:	## Do same for the consensus sequence (if '--consensus' )
				if 'cons' not in sequence_scores[HMMFname].keys():
					cons_seq = get_consensus( HMM )
					submitted = score_alignment( cons_seq , matrices, math_processes )
					cons_tally, cons_probs = sum_thread.outQ.get()
					sequence_scores[HMMFname].update( {'cons' : [cons_tally ,cons_probs]} )
				sys.stdout.write('\t{0}'.format( sequence_scores[HMMFname]['cons'][0] ) )
			sys.stdout.write('\t{0}'.format( len( seq )  ))
			sys.stdout.flush()
	##### KILL PROCESSES #####
	for i in xrange( options['-ncpus'] ):
		math_processes[i].inQ.put( 'STOP' )
	seq_send_pipe.send( 'STOP' )
	align_process.join()
	for i in xrange( options['-ncpus'] ):
		math_processes[i].join()
	calc_to_sumQ.put( 'STOP' )
	sum_thread.join()
	##### PRINT MATRIX ####
	if options['-out'] == sys.stdout:
		save = raw_input( "\nProcessed all sequences. Do you want to print the distance matrices? [y,n or filename]" )
		if re.search(r'^ye?s?$',save,re.I):
			print_HMMs_vs_seq( acc_names, sequence_scores, ali, out = sys.stdout )
		elif re.search(r'^no?$',save,re.I) or save.strip() == '':
			exit()
		else:
			sys.stdout.write( 'Saving tab-delimited data to {0}.\n'.format( save ) )
			outFile = file(save,'w')
			print_HMMs_vs_seq( acc_names, sequence_scores,ali, out= outFile)
			outFile.close()
	else:
		sys.stdout.write( '\nSaving tab-delimited data to {0}.\n'.format( options['-out'] ) )
		if os.path.exists( options['-out'] ):
			sys.stdout.write( 'Overwriting {0}.\n'.format( options['-out'] ) )
		outFile = file(options['-out'],'w')
		print_HMMs_vs_seq( acc_names, sequence_scores, ali,out = outFile )
	options['-in'].close() ## At end in case this refers to standard input, in which case raw_input won't work.
